Quantum Energy Border in Vibrational Spectra of Polynomial Molecular Potentials 
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Polynomial approximations around the global minimum of the adiabatic potential energy surface 
of several polyatomic molecules display an unphysical saddle point of comparatively small energy, 
leading to a region where the potential is negative and unbounded. This poses an energy upper limit 
for a reliable evaluation of the vibrational levels. We argue that the presence of such saddle points 
is general. The effect of tunneling through the saddle is explored by means of non-perturbative 
evaluations of vibrational spectra based on a tier by tier construction of the basis and an exact 
recursive equation for the resolvent. 

PACS numbers; 33.20.Tp, 33.15.Hp, 82.20.Kh 



I. INTRODUCTION 

The morphology of the adiabatic potential energy sur- 
face (APES), especially its low-energy minima and sad- 
dle points, is at the basis of the quantum chemistry of 
reaction paths and conformational transitions.^ The adi- 
abatic potential governs the low-energy vibrational dy- 
namics of a rigid molecule of N atoms and is a compli- 
cate function oi d — 3iV — 6 internal coordinates, such as 
bond lengths, bending or torsion angles [d = ZN — h for 
linear molecules). Its actual determination is usually a 
very demanding problemfS as the information contents of 
a function of d coordinates grows exponentially with d. 
A standard approximate parametrization is provided lo- 
cally by a Taylor expansion around the global minimum 



of these, the expansion of the vibrational Hamiltonian 
becomes: 

H = lY^huJaiP^+QD + ^^Y.'^-bcQaQbQc 
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Here, fab = V ^'^ / dqadqb\ qo is the Hessian matrix and 
yanhar collects all residual higher-order terms. For con- 
venience we take the energy zero such that V^'^{ol^) = 0. 
For a sharp minimum, the kinetic energy can be sep- 
arated into the energy of a rigid rotator, a vibrational 
term, and a mixed Coriolis term that is treated as a per- 
turbation. The simultaneous diagonalization of the Hes- 
sian matrix fab and of the vibrational kinetic quadratic 
form defines the normal-mode harmonic frequencies Ua 
and the dimensionless normal-mode coordinates Qa mea- 
sured with respect to the equilibrium geometryu^ In terms 



The frequencies uja and higher-order force constants 
^ab... of several molecules have been computed ab-initio, 
and sometimes refined by comparison with spectroscopic 
data.** For a simple molecule such as water, the force con- 
stants have been determined up to sixth order and for 
several polyatomic molecules the literature reports cal- 
culations of third and fourth-order force constantsiSiS*^ 
Hereafter, we shall refer to the (finite) polynomial poten- 
tial in Eq. Q as to the PP. 

We consider the available data for the PP of several 
molecules. In all cases we find a saddle point of compara- 
tively small energy, leading to an unphysical dissociative 
region where the potential is not bounded from below. 
We also find that, for most molecules, the potential well 
accomodates very few quantum levels up to the saddle 
energy. To our surprise, this problem seems to be ne- 
glected: its analysis is the main purpose of the present 
study. 

The standard use of the polynomial expansion ||2Jl is 
the calculation of a number of vibrational levels, usu- 
ally by means of perturbation or variational theories, or 
by numerical diagonalization. Contrary to finite-order 
perturbative calculations, a non-perturbative determina- 
tion of vibrational levels inevitably detects tunneling to 
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the unphysical region through the saddle. We show that 
the presence of this unphysical saddle is a serious issue, 
affecting radically the range of applicability of exact cal- 
culations based on the approximate PP. For this purpose 
we develop and make use of a non-perturbative technique 
(equivalent to exact diagonalization on a large basis) that 
consists in a recursive evaluation of the Green function in 
a reduced basis obtained by a "tier by tier" construction. 
The basic outcome of these calculations is that sharp vi- 
brational levels do occur only in the energy region below 
a quantum energy border, given by the sum of the energy 
of the lowest saddle (the classical border) plus a quan- 
tum correction. The latter is the zero-point energy of the 
d — 1 "transverse" modes (of positive curvature) at the 
saddle point. 

The paper is organized as follows: In Sec^lwe present 
the data for the lowest unphysical saddle point of some 
molecules, based on force constants available in the lit- 
erature, and give arguments for the general occurrence 
of a low-energy saddle in polynomial approximations of 
the molecular APES. The quantum energy border for 
the tunneling regime is introduced in Sec. IIIII The in- 
fluence of the quantum energy border on the stability 
of vibrational levels is illustrated by applying our non- 
perturbative method to a simple two-dimensional model. 
In Sec. IIVI we apply these concepts to the spectrum of 
the water molecule, and in Sec. [3 we discuss the results 
and the ensuing scenario. The computational method, 
with the main iterative formula Eq. (|A4|) , is presented in 
the Appendix. 



II. THE UNPHYSICAL SADDLE 

A unique feature of the polynomial expansion in 
Eq. (|2Jl for the special case of diatomic molecules {d—1) 
is that many even-power terms are positive (for exam- 
ple all even-power terms are positive for the Morse and 
Lennard- Jones functions). As a consequence, the trun- 
cation of the series to a positive even-power coefficient 
gives a lower-bounded potential, and thus a well-defined 
quantum problem, characterized by an infinite set of dis- 
crete levels. In addition, for the Morse potential it oc- 
curs that the second-order perturbative levels, based on 
a fourth-order PP approximation, reproduce the exact 
Morse levelsiiSi These features for d=l are unfortunately 
lost when the power expansion ^ is extended to poly- 
atomic molecules {d > 1): in all systems which we could 
obtain the anharmonic parameters for, we verified the 
occurrence of regions where the fourth-order PP is un- 
bounded below. 

Energy barriers separate different minima of a physical 
APES, corresponding to different local equilibrium con- 
figurations (isomers) of the molecule. The isomerization 
dynamics occurs mainly via quantum tunneling or ther- 
mal activation through the lowest saddle of the barrier^. 
Likewise, for the PP energy barriers separate the region 
of the physical minimum around which the expansion 




FIG. 1: A typical 1-dimensional molecular potential (Morse, 
solid) and its 5th-order polynomial approximation (dashed) 
illustrating the presence of an unbounded region separated 
from the physical confining region by a barrier topping at a 
saddle point Qs. 

is based and well grounded, from the unphysical regions 
where the potential drops to — oo. The escape to the 
unphysical region is driven by the lowest saddle which 
introduces an energy "border" that limits the range for 
(meta-stable) quantum levels allowed in the physical well 
of the PP. Figure ^ illustrates this concept in a simple 
d — 1-dimensional context, where we purposedly trun- 
cated the polynomial expansion ^ to an odd order. It 
is clear that the saddle lies in a region where the poly- 
nomial has already become a poor approximation to the 
actual APES. 

To determine the lowest saddle point of the PP of a 
polyatomic molecule, we first locate all stationary points 
in the neighborhood of Q = by solving numerically the 
equation VqF = 0, and then check that they are indeed 
saddle points, characterized by one negative and d ~ 1 
positive curvatures. Finally, for the saddle point Qs with 
lowest energy, we verify that the PP, restricted to the 
straight line through the points Q = and Qs, has a 
shape qualitatively similar to the dashed line of Fig. ^ 
i.e. that tunneling indeed occurs through a single barrier 
to a region where the potential drops to — oo. 

Table^reports the height Eg of the lowest (unphysical) 
saddle point of the PP for several polyatomic molecules, 
measured from the bottom of the potential well. The 
values of Eg of different molecules are determined by 
the characteristic anharmonicities of the molecular bonds 
and are thus fairly similar, in the 10"^ to 10^ cm"-'^ range 
for quartic PP (they are much lower for molecules char- 
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species 




£;zp(o) 




Ref. 












N2 


18619+ 


1180 


- 


12 


N2 


32834* 


1180 


- 


12 


HCl 


14919* 


1495 


- 


12 


H2O 


6846 


4717 


4803 


^ 


HOCl 


2821 


2911 


2893 


13 


H2C2O 


936 


7155 


6873 


14 


CH3OH 


50 


11398 


11276 


15 



TABLE I; Lowest saddle-point energy Es, harmonic zero- 
point energy at the minimum £'zp(0), harmonic transverse 
zero-point energy at the saddle point E^, of the 4th-order PP 
of several polyatomic molecules (*3rd or *5th-order expansion 
of Morse potential for diatomics). Energies are expressed in 
spectroscopical wavenumber units cm~^. 



actcrizcd by soft torsional modes, such as methylene 
CH3OH). Surprisingly, these saddles are low: about few 
times a typical harmonic vibrational frequency of the 
molecule. As a result, only few, if any, vibrational states 
sit below Es- This is especially true for molecules with 
a large number of atoms, where the ground-state energy 
lies substantially above the bottom of the well, by an 
amount proportional to d, that easily compares or ex- 
ceeds the saddle energy Eg. 

The occurrence of a saddle leading to an unphysical 
region is by no means specific of the PP of the molecules 
listed in Table we argue that this is a general feature 
to be found in the PP of most polyatomic molecules. 
Indeed, also for an even-power truncation, the PP can 
easily drop to —00 in some direction in Q spacc"^^. The 
precise value of the 4th-order parameters ^abcd (includ- 
ing their sign) is determined by the local properties of 
the physical APES at its minimum Q = 0, not by any re- 
quirement of confining behavior at large distance: the PP 
of a real molecule easily contains negative semi-diagonal 
terms ^aabb and sizable mixed terms ^abcd, which in turn 
produce regions where the PP drops to —00. In prac- 
tice, the same argument prevents confining behavior also 
of 6th and higher even-order terms, and an even-power 
truncated PP does not behave any better, away from 
the physical minimum, than an odd-power truncated 
PP. Therefore, we consider it extremely unlikely (though 
technically possible) that a real polyatomic molecular po- 
tential may ever be found whose polynomial expansion 
at the minimum (truncated at any order > 2) is lower- 
bounded everywhere. 



III. THE QUANTUM ENERGY BORDER 

As the PP has no lower bound, the associated 
Schrodinger problem is ill-defined. However, resonant 
states with complex energy values Ea — iTa usually exist 
in the well. The imaginary part is proportional to the 



inverse lifetime of the state, and it is easily seen to be 
also proportional to the tunneling current through the 
barrier. Pa can be evaluated by semiclassical methods, 
like WKB or instanton approacbiSiiLi^, or by complex 
dilatation^^'^*^'^*. Deep in the well, Ta is in general ex- 
ponentially small and proportional to the Gamow factor 
g-2S/s^ where S is the imaginary-time action along the 
most probable tunneling path inside the barrier. Reso- 
nances appear as sharp peaks in the spectral density, at 
energies Ea extremely close to the eigenenergies of the 
Schrodinger equation restricted to the well. As energy 
increases toward the threshold value, tunneling evolves 
into a "leaking" regime, characterized by the appearence 
of resonances whose imaginary part Fq is comparable to 
the real parlsW^ relics of further excited states in the well, 
strongly hybridized to the continuum. In one dimension, 
the threshold energy coincides with the saddle one Eg. 

For d > 1, an effect absent in the simple one- 
dimensional picture of Fig. ^ is to be considered: tunnel- 
ing through the barrier at the saddle point is hindered 
by the "transverse" motion of the degrees of freedom per- 
pendicular to the one crossing the barrier. These perpen- 
dicular degrees of freedom are associated to a minimum 
energy E^p due to Heisenberg's uncertainty, that adds 
to the saddle height to determine the quantum energy 
border between the tunneling and the leaking regimes 



E, 



qb 



Eg + Ei 
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As the study of tunneling problem: 
proximate E:^^ by its harmonic expression 



(3) 



suggests, we ap- 



E. 



zp 



1 



(4) 



in terms of the d — 1 real harmonic frequencies oj[ at 
the saddle point {uji is the imaginary frequency associ- 
ated to the tunneling direction). Since, for most poly- 
atomic molecules, i?^ is close to the zero-point energy 
of the ground state _Ezp(0) (see Table QJ, the raising of 
the energy border due to recovers a spectral range 
(i?zp(0) ^ E < Eqb) where quasistationary vibrational 
levels are to be found, of extension comparable to the 
classical region for bounded motion {0 < E < Eg). 

To illustrate the quantum energy border, we propose 
and discuss the following two-dimensional model Hamil- 
tonian 



H 

U{x,y) 



^ + U{x, y) 



(5) 




where the vibrational coordinates x and y are dimension- 
less, like Qa in Eq. Q. This potential is simple enough to 
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FIG. 2: The polynomial potential U(x,y) of Eq. We 
take hwx as the basic energy units in this model. The other 
parameters are: ujy = 1.7 ujx, Es = 7 tiuox, with (a) lus = uix, (b) 
ujs — 2.7 UJ2:, and (c) u}s~?!i-^x- 



allow for an analytic control of its shape through the ad- 
justable parameters uj^, t^Jy, lOs and Eg. Figure [21 depicts 
U{x, y) for a few choices of these parameters. The origin 
is a local minimum, with harmonic frequencies ujx and 

LUy. A saddle point {xs, Us) — ( \/(yEs/fkJx, ) , of height 



Es above the minimum, separates the stable region from 
the region of large positive x, where the potential has no 
lower bound. The quadratic expansion of the potential 
near the saddle point yields an imaginary frequency for 
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FIG. 3: The spectrum of the model of Eq. parameters 
as in Fig.|5Jb), initial state \vx =5,Vy =0) , and e — 0.005 hujx- 
Bottom to top: N=10 (396 states), iV=50 (7956 states), 
A^=100 (30906 states). The sequence of converged peaks rep- 
resents the Vx progression. Dense satellite lines (not con- 
verging with TV) slightly below and above the quantum en- 
ergy border (dashed line) Eq^ indicate tunneling and leakage 
broadening. 



the X motion uj'^ = ioj^, and a harmonic one ui'y = ujs 
for the surviving transverse mode. The y zero-point en- 
ergy E^p — tiuJs/2 extends the tunneling regime up to the 
quantum energy border Eqh = Es + hujs/2. Sharp vibra- 
tional levels should therefore be found in the approximate 
range {fujx + fiujy)/2 ^ E < Eg + fujJs/2. 

These expectations have been checked numerically by 
computing the eigenvalue spectrum. A standard way to 
evaluate the eigenvalues of an operator H is to compute 
the Green function 



G{E + ie) = (vo|(S + i£-i7)-i|vo) 



(6) 



for real E and a small positive imaginary part e. The 
eigenvalues Ea of the exact eigenstates \a) of the Hamil- 
tonian appear as poles in the complex plane slightly be- 
low the real axis. The reference state |vo) represents the 
initial excitation; for our model, Eq. (|SJ, we take for |vo) 
an eigenstate \vx, Vy) of the harmonic part of the Hamil- 
tonian. The imaginary part of —G{E + ie)/'iT yields the 
lineshape function 



m 



E 



(a|vo) 



{E - EaY + £2 



(7) 



Eigenvalues show up as peaks of I{E), with heights pro- 
portional to the squared overlap of the exact eigenstates 
to the initial excitation. The parameter e introduces an 
artificial Lorentzian broadening of the lines Ea- 
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FIG. 4: A blowup of the spectrum of Fig. |3 in the region 
near -Bqb, computed on bases with A'' = 10, 11, ... 20 tiers. 
The strongly basis-dependent satellites build up an envelope 
representative of tunneling and leakage broadening (dashed 
line). The envelope width of the state close to 7 hu! is much 
smaller than the resolution e = 0.005 huix chosen for this cal- 
culation. The width of the next peak, at 7.8 hw, is estimated 
r ~ O.Oihuj; for the structures above i?qb, broadening be- 
comes even larger, comparable to the inter-level separation. 



The method for evaluating the Green function is non- 
perturbative, recursive, and it is applied in two steps. 
First, one constructs a basis of harmonic states grouped 
into families Ti (tiers) , i = . . . N, adapted to the spe- 
cific PP and to the choice of the initial "bright" state 
|vo) (which constitutes tier J|j) 23i24,25 _ Next, the func- 
tion G{E + ie) is evaluated through the exact recursive 
relation Eq. l|A4p . as detailed in the Appendix. Conver- 
gence of the spectrum against the basis size is tested by 
changing the single parameter N , the number of tiers. 
As the basis is finite, the spectrum consists of a finite 
number of real eigenvalues. Nevertheless, tunneling and 
leakage linewidth broadening shows in the comparison of 
spectra I{E) of different tier number N . As we shall see, 
eigenvalues well below the quantum energy border Eq\^ 
correspond to stable sharp peaks. Near the border and in 
the leaking regime, evaluations with different N produce 
a number of weaker unphysical poles (not converging for 
increasing N) that decorate the resonant peaks. 

Figure|31reports the spectra computed for three largely 
different bases for the potential of Fig.l^^b) {uJs = 2.7 lOx). 
Stable sharp levels of the Vx progression occur until 
slightly below the quantum border Eq^, — 8.35 hujx. 
Starting close to £^qb, dense non-converging satellites 
spring up, indicating a tunneling broadening exceeding 
the arbitrarily chosen spectral resolution e — 0.005 hujx. 
The proliferation of spurious peaks, which marks the on- 
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FIG. 5: The relative difference between the spectra of Fig. El 
computed on two very different bases (A'^ = 50 and A'' = 100) 
as a function of energy for the three potentials of Fig. |5| The 
narrowing of the saddle valley - with associated increase in 
iSqb = 7.5huJx, 8.35 /ioJi:, and llhujx (vertical dashed lines) - 
determines an enlargement of the region of convergence. 



set of the leaking regime, is more manifest if one super- 
poses several evaluations of I{E) for different values of N . 
Figure 0] is a magnification of Fig. |31 in the crossover re- 
gion near i?qb. The fitting envelope exhibits sharp peaks 
in the tunneling regime and substantial broadening in the 
leaking regime. 

To better illustrate the role of = hujs/2, Fig.Oplots 
the relative difference between the spectra computed on 
two different bases {N = 50 and N = 100) for the three 
potentials of Fig. [3 In the energy region of weak tun- 
neling, this relative difference, associated to the weak 
satellite peaks representative of broadening, is extremely 
small, and it grows roughly exponentially with the en- 
ergy distance to £^qb- Above i?qb, the relative difference 
is close to unity, indicating strong leakage. As expected, 
the region of weak tunneling enlarges with fiwg, following 
the increasing Eqi,. 

We recall that the spectra of Fig. O and [S] correspond 
to an initial excitation given by the fourth co^ overtone 
\vx — 5,Vy = 0), a state with a substantial extension in 
the direction of the saddle. If a "transverse" initial exci- 
tation such as \vx = 0,Vy = A) is considered instead (see 
the spectrum in Fig.|^, tunneling is strongly suppressed, 
since the overlap of the initial excitation with the fast- 
leaking states is very small. Indeed, the main spectral 
structures are affected by much smaller broadening, even 
above the saddle, in comparison to the spectra of Figs. El 
and 1^1 This comparison exhibits a general feature of 
multidimensional PP (o?> 1), namely the coexistence of 
weakly and strongly leaking states in the same energy 
range (above -Eqb). For large d, the latter represent a 
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FIG. 6: Convergence of the spectrum as in Fig. |3] but for 
initial excitation \vx~0,Vy — 4) perpendicular to the direction 
of the saddle. Below i?qb (dashed line), and even above it, the 
main features are affected by much smaller broadening (few 
weak non-converging satellites, a substantially smaller relative 
difference in the lower panel) than with initial excitation \vx = 
5,vy^0) of Figs.OEl 

minor fraction of all the states because of the small frac- 
tion of leaking directions over the whole solid angle in 
configuration space. 



IV. A REAL MOLECULE: H2O 

As an illustration of the applicability of the tiers 
Green-function method to a real molecule and of the 
difficulties brought in by the energy border of the PP, 
we compute the vibrational spectrum of water. We em- 
ploy the ab-initio quartic PP parameters listed in Table 
6 of Ref. 1^ The two real frequencies huj'2 = 4758 cm^^ 
and Swg = 4849 cm^^ at the lowest saddle point of H2O 
{Es = 6846 cm~^) are found so much larger than the 
three curvatures at the minimum, that E^p > E^p{0) 
(Table nj: this produces a rather narrow saddle, which 
pushes the quantum border up to Eqi, = 11649cm~^. 

We take the harmonic fundamental excitation of the 
uji (symmetric OH stretching) mode as initial state, and 
obtain the spectrum in Fig. \7\ The peaks, represent- 
ing exact vibrational levels, are assigned to the harmonic 
quantum numbers of the closest state resulting from stan- 
dard second-order perturbation theory. Convergence is 
studied by increasing the tier number from TV = 5 (not 
shown in figure) to = 10 and 15: it is very good al- 
ready using N — 5 tiers, thus showing the effectiveness 




4000 6000 8000 10000 12000 14000 
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FIG. 7: The spectrum of H2O for initial excitation = 
1,V2 =0,113 = 0), computed with A''=10 (5570 states) and 
N=15 (16811 states), broadening e—i cm~^. Below the quan- 
tum energy border (dashed line) convergence is satisfactory 
with relative error < 1%. Like in Fig. |S] the main features 
seem to converge also above -Bqb, but the appearance of new 
structures for larger A'^ indicates non-negligible leakage. 



of the tiers Green-function method applied to a realis- 
tic problem. Like in the 2-dimensional model of Sec. IIIII 
hardly any tunneling is observed below Eq^. Across -Eqb, 
the appearence of the weak non-converging satellites con- 
firms that the resonances in this region are affected by 
significant leakage to the continuum. Indeed, since the 
lowest saddle lies in a direction involving mainly mode 1, 
the chosen initial state |1, 0, 0) encourages tunneling, like 
in the example of Fig. O Nonetheless, the stabililty of 
several features even above Eq^ indicates that a number 
of fairly long-lived physical states are so localized in the 
well that their overlap to the outside continuum is rela- 
tively small. We are currently investigating a method to 
obtain an explicit evaluation of the tunneling widths Pq 
of the resonance states, by means of a complex coordi- 
nates transformation inspired by Ref. Il9l 



V. DISCUSSION 

We show that, in general, polynomial approximations 
of molecular potentials display unphysical saddles that 
lead to regions where the potential is not lower-bounded. 
The height of the saddle and the zero-point energy of 
transverse modes both determine a quantum energy bor- 
der. 

In the spectral region below the quantum border, tun- 
neling is exponentially small, the PP is reliable and stan- 
dard perturbative treatment of the anharmonic interac- 
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tions usually provides a good approximation to the nar- 
row level positions. Perturbative calculations are often 
extended to higher energies, based on the general belief 
that, like in the d—1 Morse case, the second-order ap- 
proximation compensates the wrong behavior of the PP 
away from the minimum and reproduces the levels of a 
physical APESifi. 

Our non-perturbative calculations show that, as energy 
is raised above the quantum border, resonances leave the 
tunneling regime and couple more and more strongly to 
the continuum, practically washing out all spectral de- 
tail. Such highly excited states are of scarce physical 
relevance anyway, since the extension of the associated 
wavefunction explores regions where the PP becomes a 
very poor approximation of the true APES. 

The present results cast a shadow on the traditional 
use of the polynomial approximation of the physical 
APES for the calculation of highly excited vibrational 
or roto-vibrational spectra of polyatomic molecules. For 
example, in most molecules, IVR (intramolecular vi- 
brational relaxation) spectra involve energy levels much 
above the saddle. Hence, IVR applications of exact 
numerical methods, such as the Lanczos or Davidson 
algorithma^SiSLSSiSS or the Green-function method pro- 
posed here, are bound to face the problem of tunneling. 
The consequent broadening of the resonant states poses 
an intrinsic limit to the spectral accuracy which can be 
obtained. A real progress may be achieved by the use 
of smarter parametrizations of the APES (e.g. based on 
Morse coordinates), providing the correct large-Q behav- 



APPENDIX A: NON-PERTURBATIVE 
EVALUATION OF THE GREEN FUNCTION. 



We present a general procedure, inspired to Ref. |30 
and there indicated as "tiers method" , for the non- 
perturbative evaluation of the eigenvalues and spectral 
weights of a Hamiltonian decomposed as H = Hq + V. 
It is assumed that the eigenvectors of Hq are well char- 
acterized and that the matrix elements of V are easy 
to compute and link any eigenstate of Hq with a finite 
(not too large) number of other eigenstates. For defi- 
niteness, we consider the problem at hand of the vibra- 
tional levels in a polynomial potential The term 
Hq = ^ HLUi{alai + 1/2) describes d independent oscilla- 
tors, with eigenvectors |v) ~ \vi, V2, ■ ■ ■ , Vd) that specify 
the occupation numbers of the d oscillators. V is the 
anharmonic part of the potential, yanhar^ 

The first step of the method is to partition the unper- 
turbed eigenvectors in families (tiers), Tq,Ti,... of de- 
creasing perturbative relevance. The construction is such 
that the matrix representation of H is block-tridiagonal 
in the tiers. Symbolically we write the blocks as Ha = 
{T,\H\Ti) and = {T,\H\T,+i). 

Next, we provide an exact iteration scheme for the eval- 
uation of the matrix elements of the resolvent G{z) = 



[zl — H)~^ in the basis vectors |vo,q) in Tier 0. The 
poles give the eigenvalues Ea of iJ, and the residues are 
the spectral weights: 



\Ea){E 
Z-Ea 



(Al) 



The actual evaluation requires approximations such as 
truncation of the tier generation to some order, and 
possibly restriction of tier population. Anyhow, as the 
method is iterative and in analogy with procedures that 
approximate a self-energy, the resulting matrix elements 
of G^"-' (z) are highly non-perturbative. 

Tier construction 
Depending on the problem under investigation, a set of 
to unperturbed states |vo,q) {a = l,...,to) is selected 
to form the initial tier Tq. In the computations of this 
paper, Tq contains a single state, the bright state. The 
action of V on Tq gives new vectors; the basis states that 
have non-zero overlap with one of them, and are not in 
To, are collected in Ti. We label them as |vi,q) [a — 
1 . . . ti). For a finite set To, and an interaction V which 
is a polynomial in the raising and lowering operators, tier 
Ti and subsequent ones are finite. 

Next we consider the set VTi^ and expand it in the 
eigenvectors already in To and Ti , plus new ones that are 
collected in tier T2. The process is iterated to produce 
further tiers T3, T4, ... Up to this point the method is sim- 
ply a smart algorithm to generate systematically a good 
approximate basis for a quantum problem where some 
non-interacting part Hq of the Hamiltonian can be sin- 
gled out. Indeed, such a basis has been employed success- 
fully in a different context, in conjunction with standard 
(Lanczos) techniques^ However, the hierarchical basis 
structure and the corresponding block-tridiagonal form 
of the Hamiltonian, suggest a natural iterative method 
to construct the spectrum. 

Evaluation of the Resolvent 
The iterative method is based on the following formula 
for the inversion of partitioned matrices, with square di- 
agonal blocks (we omit unneeded terms): 



M 



M' 



Mil Mi2\ 
^M21 M22y 

'[Mli-Mi2(M22)-lAf2l] 



(A2) 



We apply this formula by identifying M with the matrix 
representations oi zl ~ H = {zl — Hq) — V and with 
G{z). The blocks result from the separation of the basis 
into the set Tq and the ordered set Ti U T2 U . . . In this 
tier basis, off-diagonal matrix elements of M only arise 
from the action of V . Thus Mn ~ zIq — Hqq {Iq is the 
unit matrix of size to and Hqq — (To|-ff |To)). M22 is the 
matrix {zl — H) expanded in the remaining tiers. M12 = 
M|i is a rectangular matrix of size to x (ti -I- 12 + . . .); 
by the tier construction, non-zero matrix elements of the 
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potential are only in the leftmost submatrix of size toxti, 
that identifies with the matrix — Vqi = — (To|y|ri). 

The aim of the calculation is to evaluate the block 
(Af~i)ii EE G'(°)(z), Eq. (HU. The inversion formula 
(|X2|) provides 

G(°)(z) = [Afn-Afi2(M22)-^Af2i]-i 

= [zio - Hoo - VoiG^'Hz)Vwr' (A3) 

We defined thehxti matrix G'(i)(z) = (ri|(M22)~^|Ti). 
To evaluate it we use Eq. (IA2|) again, with the blocks now 
resulting from the separation of the basis into the set Ti 
and the set T2Ur3U. . .. Now the block Mn is {zIi-Hn) 
and iV/22 is the matrix {zl — H) expanded in the basis 
r2 U T3 U . . . The matrix G^^^ (z) coincides with the block 
(M~i)ii. The formula IIA2II provides: 

G(i)(z) = [zh - ffii - V^i2G(2)(z)V^2i]-^ 

where G(2)(z) = (r2|(Af22)"^|T2). By iterating the same 
inversion formula IIA2p we obtain a chain of relations of 
the type 

G(^-i)(z) = [zh-i-Hk-i,k-i-Vk-i.kG^''\z)Vk.k-ir^ 

(A4) 

In practice, the (in principle infinite) chain is truncated 
by approximately taking G^^\z) « {zIn — Hffj^)~^. 
This assumption is in order if the coupling of Tjv to 
the subsequent tier is negligible. Starting from G^^'(z), 
one iterates ljA4|l back to the sought for matrix G'-'^'(z). 
This procedure is a matrix generalization of the contin- 
ued fraction expansion for the inversion of tridiagonal 
matrices 1^ 
Discussion 

The method just outlined generates an effective basis and 
solves the quantum mechanical problem in the reduced 
Hilbert space exactly. The basis can be enlarged system- 
atically by increasing the number of tiers TV, to achieve 



convergence. When applied to a lower-bounded PP, this 
method provides an efficient and rapidly converging nu- 
merical scheme to determine its exact eigenenergies. For 
the lower-unbound PP of a polyatomic molecule, this 
method provides good estimates of the position of the 
quasi-stationary states, including a rigorous treatment of 
all anharmonic resonances. The recursive calculation of 
the Green's function (jA4p has several advantages with re- 
spect to the more traditional Lanczos methodiSLSLSiM 
(i) it provides equal accuracy through the whole spec- 
trum, while Lanczos method is more accurate close to the 
endpoints; (ii) it splits the Hilbert space into subspaces 
To, ...Tat to treat one at a time; (iii) once the chain of 
matrices is set up, each frequency requires an indepen- 
dent calculation, which makes this method suitable for 
parallel calculations. Its main disadvantage is the rapid 
growth of the tier size ti, for systems with many degrees 
of freedom. To fit the available CPU/memory limits, it is 
possible to cutoff the tier growth to some maximum size 
^max, as described in Ref. [sj- In general, the recursive 
method may become very costly in CPU time, since the 
evaluation of G^"-* {E + ie) requires N inversions for each 
sample frequency E, each inversion costing a time pro- 
portional to i^ax- the Lanczos method, a single chain 
of A'^Lanczos « 10^ iterations, each costing of the order of 
the Hilbert space size ^ N ■ imax, generates the whole 
spectrum. 

The C++ code for computing the tier basis and the spec- 
trum based on the Green-function recursive inversion for- 
mula ljA4p is available in Ref. 1351 
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